/*
 * Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.
 *
 * Please refer to the NVIDIA end user license agreement (EULA) associated
 * with this source code for terms and conditions that govern your use of
 * this software. Any use, reproduction, disclosure, or distribution of
 * this software and related documentation outside the terms of the EULA
 * is strictly prohibited.
 *
 */

// Simple 3D volume renderer

#ifndef _VOLUMERENDER_KERNEL_CU_
#define _VOLUMERENDER_KERNEL_CU_

#include <cutil_inline.h>    // includes cuda.h and cuda_runtime_api.h
#include <cutil_math.h>

#include "volumeRender_kernel_common.cuh"


/**
 *	d_render(uint *d_output, uint imageW, uint imageH,
 *             float density, float brightness,
 *             float transferOffset, float transferScale, 
 *             float kgc,float kgs, float kge)
 *	 
 *	 Render it.
 */
__global__ void
d_render_tone(uint *d_output, uint imageW, uint imageH,
         float density, float brightness,
         float transferOffset, float transferScale, 
		 float kgc,float kgs, float kge)
{
    const int maxSteps = 500;
    const float tstep = 0.01f;
    const float opacityThreshold = 0.95f;
    const float3 boxMin = make_float3(-1.0f, -1.0f, -1.0f);
    const float3 boxMax = make_float3(1.0f, 1.0f, 1.0f);

	uint x = blockIdx.x*blockDim.x + threadIdx.x;
    uint y = blockIdx.y*blockDim.y + threadIdx.y;
    if ((x >= imageW) || (y >= imageH)) return;

    float u = (x / (float) imageW)*2.0f-1.0f;
    float v = (y / (float) imageH)*2.0f-1.0f;

    // calculate eye ray in world space
    Ray eyeRay;
    eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
    eyeRay.d = normalize(make_float3(u, v, -2.0f));
    eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

    // find intersection with box
	float tnear, tfar;
	int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
    if (!hit) return;
	if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

    float4 sum = make_float4(0.0f);
    float t = tnear;
    float3 pos = eyeRay.o + eyeRay.d*tnear;
    float3 step = eyeRay.d*tstep;

    /*The variables needed for tone shading*/
    float ktd=0.6f;
    float kta=1.0;
    float ktb=0.3f;
    float kty=0.3f;
    
    float4 cw;
    cw.x=kty;
    cw.y=kty;
    cw.z=0.0f;
    cw.w=1.0f;

    float4 cc;
    cc.x=0.0f;
    cc.y=0.0f;
    cc.z=ktb;
    cc.w=1.0f;

    float4 It;
    float4 Io;
    float3 light_orig;
    light_orig.x=1;
    light_orig.y=1;
    light_orig.z=1;
    float3 light;
    float light_mag=0;
 

   for(int i=0; i<maxSteps; i++) {
	    
        float posx=  pos.x*0.5f+0.5f;
	    float posy=  pos.y*0.5f+0.5f;
	    float posz= pos.z*0.5f+0.5f;
        
	    float sample = tex3D(tex, posx, posy, posz);
	    float gradMag = tex3D(gradMagTex, posx, posy, posz);
        float4 gradvec = tex3D(gradientsTex, posx,posy,posz);
	    float4 col = tex1D(transferTex, (sample-transferOffset)*transferScale);

        //col.w *= density;
       
		/*Tone Shading*/
		light.x=abs(pos.x-light_orig.x);
		light.y=abs(pos.y-light_orig.y);
		light.z=abs(pos.z-light_orig.z);

		light_mag=sqrtf(light.x*light.x+light.y*light.y+light.z+light.z);
		light.x/=light_mag;
		light.y/=light_mag;
		light.z/=light_mag;
		
		float lgrad=abs(gradvec.x*light.x) + (gradvec.y*light.y) + (gradvec.z*light.z);

	   	It.x= ((1.0+lgrad)/2)*cw.x + (1- (1.0+lgrad)/2)*cc.x;
	   	It.y= ((1.0+lgrad)/2)*cw.y + (1- (1.0+lgrad)/2)*cc.y;
	   	It.z= ((1.0+lgrad)/2)*cw.z + (1- (1.0+lgrad)/2)*cc.z;
		if(lgrad>0)
		{
			Io.x= ktd*1*lgrad;
			Io.y= ktd*1*lgrad;
			Io.z= ktd*1*lgrad;
		}
		else
		{
			Io.x=0;
			Io.y=0;
			Io.z=0;
		}
		col.x= It.x+ktd*Io.x;
		col.y= It.y+ktd*Io.y;
		col.z= It.z+ktd*Io.z;
		//col.x=gradvec.x;
		//col.y=gradvec.y;
		//col.z=gradvec.z;


        //pre-multiply alpha
        col.x *= col.w;
        col.y *= col.w;
        col.z *= col.w;
        
        sum = sum + col*(1.0f - sum.w);  // "over" operator for front-to-back blending

        if (sum.w > opacityThreshold)
            break;

            t += tstep;
            if (t > tfar) break;

            pos += step;
    }
    sum *= brightness;

    // write output color
    d_output[y*imageW + x] = rgbaFloatToInt(sum);
}

//


#endif // #ifndef _VOLUMERENDER_KERNEL_CU_
